/* 9/11/2013 */
/* basic models */

#delimit;
set more off; 
/*log using 9_11_13.log, t replace; */

use coaldata_BJPS;

/* model 1 */

logit govtparty incumb_parties largest_party seatshare seatshare_sq 
prevpm distance_median_weighted if numparties==1, cluster(formnum);

matrix coeffs_l = e(b);
matrix covmat_l = e(V); 

/* model 2 */

clogit govt_coal incumb_parties largest_party seatshare seatshare_sq prevpm distance_median_weighted, group(formnum);

/* model 3 */

clogit govt_coal incumb_parties largest_party seatshare seatshare_sq prevpm distance_median_weighted incumb_govt, group(formnum);

/* model 4 */
/*
clogit govt_coal incumb_parties incumb_govt largest_party seatshare 
seatshare_sq prevpm maxid_diff minor minwin numparties, group(formnum);
*/
clogit govt_coal incumb_parties largest_party seatshare seatshare_sq prevpm incumb_govt 
distance_median_weighted minor minwin numparties, group(formnum);

predict cp, p;

/* McFadden train test */

sort formnum;

by formnum: egen incumb_partiesS = sum(incumb_parties * cp);
by formnum: egen incumb_govtS = sum(incumb_govt * cp);
by formnum: egen largest_partyS = sum(largest_party * cp);
by formnum: egen seatshareS = sum(seatshare * cp);
by formnum: egen seatshare_sqS = sum(seatshare_sq * cp);
by formnum: egen prevpmS = sum(prevpm * cp);
by formnum: egen distance_median_weightedS = sum(distance_median_weighted * cp);
by formnum: egen minorS = sum(minor * cp);
by formnum: egen minwinS = sum(minwin * cp);
by formnum: egen numpartiesS = sum(numparties * cp);

generate incumb_partiesZ = 0.5*(incumb_parties - incumb_partiesS)^2;
generate incumb_govtZ = 0.5*(incumb_govt - incumb_govtS)^2;
generate largest_partyZ = 0.5*(largest_party - largest_partyS)^2;
generate seatshareZ = 0.5*(seatshare - seatshareS)^2;
generate seatshare_sqZ = 0.5*(seatshare_sq - seatshare_sqS)^2;
generate prevpmZ = 0.5*(prevpm - prevpmS)^2;
generate distance_median_weightedZ = 0.5*(distance_median_weighted - distance_median_weightedS)^2;
generate minorZ = 0.5*(minor - minorS)^2;
generate minwinZ = 0.5*(minwin - minwinS)^2;
generate numpartiesZ = 0.5*(numparties - numpartiesS)^2;

/* McFadden-Train LM test */

clogit govt_coal incumb_parties incumb_govt largest_party seatshare 
seatshare_sq prevpm distance_median_weighted minor minwin numparties 
incumb_partiesZ incumb_govtZ largest_partyZ seatshareZ 
seatshare_sqZ prevpmZ distance_median_weightedZ minorZ minwinZ numpartiesZ, group(formnum);


mixlogit govt_coal largest_party seatshare seatshare_sq distance_median_weighted prevpm minor, 
rand(incumb_parties incumb_govt minwin numparties) 
nrep(200) group(formnum);

matrix coeffs = e(b);
matrix covmat = e(V);

keep if formnum==155; /* Iceland 1995 */
set obs 1000;

/* random draws of coefficients */

drawnorm b1-b14, means(coeffs) cov(covmat) double;


generate incumb_parties2 = incumb_parties;
generate incumb_govt2 = incumb_govt;


forvalues k = 1/6 {;
  generate baseprob`k' = .;
  generate cfprob`k' = .;
  generate diff`k' = .;
};

/* seatshares: 
party2 = 4.76
party5 = 23.81
party6 = 39.68
party1 = 14.29 
party4 = 11.11
party3 = 6.35 

6 and 4 were prev govt, switch that to 6 and 5
*/

nois _dots 0, title() reps(1000);
quietly forvalues j = 1/1000 {;

  mkmat b1-b14 if _n==`j', matrix(bdraw);
  mixlpred2 cabprob_base; 

  qui replace incumb_parties = 0 if partydum4==1;
  qui replace incumb_parties = 1 if partydum5==1;
  qui replace incumb_govt = 0;
  qui replace incumb_govt = 1 if partydum1==0 & partydum2==0 & partydum3==0 & partydum4==0 & partydum5==1 & partydum6==1;

  mixlpred2 cabprob_cf; 
  
  qui replace incumb_parties = incumb_parties2;
  qui replace incumb_govt = incumb_govt2;
  
 
  
  forvalues n = 1/6 {;
  egen t_baseprob`n' = total(cabprob_base) if partydum`n' == 1;
  egen baseprob_temp`n' = max(t_baseprob`n');
  egen t_cfprob`n' = total(cabprob_cf) if partydum`n' == 1;
  egen cfprob_temp`n' = max(t_cfprob`n');
  
  replace baseprob`n' = baseprob_temp`n' if _n==`j';
  replace cfprob`n' = cfprob_temp`n' if _n==`j';
  replace diff`n' = cfprob`n' - baseprob`n' if _n==`j';
  
  drop t_baseprob`n' baseprob_temp`n' t_cfprob`n' cfprob_temp`n';
  };
  
drop cabprob_base cabprob_cf;

nois _dots `j' 0;
  
};





/* logit counterfactual */

drop incumb_parties2;

drawnorm bl1-bl7, means(coeffs_l) cov(covmat_l) double;

/* baseline and counterfactual profiles for each party */
generate incumb_parties1 = 0;
generate largest_party1 = 0;
generate seatshare1 = 14.28571;
generate seatshare_sq1 = 2.040815;
generate prevpm1 = 0;
generate distance_median_weighted1 = .04325;

generate incumb_parties2 = 0;
generate largest_party2 = 0;
generate seatshare2 = 4.761905;
generate seatshare_sq2 = .2267574;
generate prevpm2 = 0;
generate distance_median_weighted2 = .00655;

generate incumb_parties3 = 0;
generate largest_party3 = 0;
generate seatshare3 = 6.349206;
generate seatshare_sq3 = .4031242;
generate prevpm3 = 0;
generate distance_median_weighted3 = .08885;

generate incumb_parties4 = 1;
generate largest_party4 = 0;
generate seatshare4 = 11.11111;
generate seatshare_sq4 = 1.234568;
generate prevpm4 = 0;
generate distance_median_weighted4 = .00655;

generate incumb_parties5 = 0;
generate largest_party5 = 0;
generate seatshare5 = 23.80952;
generate seatshare_sq5 = 5.668933;
generate prevpm5 = 0;
generate distance_median_weighted5 = .00795;

generate incumb_parties6 = 1;
generate largest_party6 = 1;
generate seatshare6 = 39.68254;
generate seatshare_sq6 = 15.74704;
generate prevpm6 =  1;
generate distance_median_weighted6 = .30575;

generate incumb_parties1c = 0;
generate incumb_parties2c = 0;
generate incumb_parties3c = 0;
generate incumb_parties4c = 0;
generate incumb_parties5c = 1;
generate incumb_parties6c = 1;

forvalues v = 1/6 {;
generate baseprob_l`v' = exp(bl7 + incumb_parties`v'*bl1 +  largest_party`v'*bl2 + 
seatshare`v'*bl3 + seatshare_sq`v'*bl4 + prevpm`v'*bl5 + 
distance_median_weighted`v'*bl6)/
(1+exp(bl7 + incumb_parties`v'*bl1 +  largest_party`v'*bl2 + 
seatshare`v'*bl3 + seatshare_sq`v'*bl4 + prevpm`v'*bl5 + 
distance_median_weighted`v'*bl6));

generate cfprob_l`v' = exp(bl7 + incumb_parties`v'c*bl1 +  largest_party`v'*bl2 + 
seatshare`v'*bl3 + seatshare_sq`v'*bl4 + prevpm`v'*bl5 + 
distance_median_weighted`v'*bl6)/
(1+exp(bl7 + incumb_parties`v'c*bl1 +  largest_party`v'*bl2 + 
seatshare`v'*bl3 + seatshare_sq`v'*bl4 + prevpm`v'*bl5 + 
distance_median_weighted`v'*bl6));

generate diff_l`v' = cfprob_l`v' - baseprob_l`v';

};

forvalues m = 1/6 {;
summarize diff`m';
generate sddiff`m' = r(sd);
summarize diff_l`m';
generate sddiff_l`m' = r(sd);

generate diffdiff`m' = diff`m' - diff_l`m';
generate sddiffdiff`m' = sqrt((sddiff`m'^2) + (sddiff_l`m'^2));
generate zdiffdiff`m' = abs(diffdiff`m'/sddiffdiff`m');
};

/* Party Numbers:
1: People's Alliance
2: Women's List 
3: National Awakening
4: Social Democratic
5: Progressive 
6: Independence
*/

/* mixed logit counterfactual probabilities */
summarize 
baseprob1 cfprob1 diff1 
baseprob2 cfprob2 diff2 
baseprob3 cfprob3 diff3 
baseprob4 cfprob4 diff4 
baseprob5 cfprob5 diff5 
baseprob6 cfprob6 diff6;

/* binary logit counterfactual probabilities */
summarize 
baseprob_l1 cfprob_l1 diff_l1 
baseprob_l2 cfprob_l2 diff_l2 
baseprob_l3 cfprob_l3 diff_l3 
baseprob_l4 cfprob_l4 diff_l4 
baseprob_l5 cfprob_l5 diff_l5 
baseprob_l6 cfprob_l6 diff_l6;

/* test for differences in differences */
/* Significant if abs(z) > 1.96 */

summarize 
zdiffdiff1 zdiffdiff2 zdiffdiff3 
zdiffdiff4 zdiffdiff5 zdiffdiff6;

